rm(list = ls())

library(causalsens)
library(DirectEffects)
library(haven)
library(foreign)
library(rio)
library(tidyverse)
library(sjlabelled)
library(memisc)
library(expss)
library(dplyr)
library(tidyverse)
library(magrittr)
library(sandwich)
library(miceadds)
library(stargazer)
library(lfe)

cw <- import("/Users/christopherblair/Desktop/JOP Replication/ams175_R.dta")

##### SEQUENTIAL G -- MAIN

#ESTIMATE CHAMBER COEFFICIENT

chamber1_prepost <- lm(use ~ numchamber+ nummask+
              age2+age3+age4+age5+age6+age7+age8+age9+age10+age11+age12+school2+school3+school4+school5+school6+rankgrade2+rankgrade3+rankgrade4+rankgrade5+rankgrade6+monthsoversea2+monthsoversea3+monthsoversea4+monthsoversea5+monthsoversea6+monthsoversea7+monthsoversea8+monthsoversea9+monthsoversea10+base1+base2+base3+base4+
              ally_relations+orientation_officers+infocenter+orient_meet+war_interest+honolulu_contact,
              data=cw)
summary(chamber1_prepost, cluster="base")

chamber1_g <- lm(I(use -
              coef(chamber1_prepost)['ally_relations']*ally_relations +
              coef(chamber1_prepost)['orientation_officers']*orientation_officers +
              coef(chamber1_prepost)['orient_meet']*orient_meet +
              coef(chamber1_prepost)['war_interest']*war_interest) ~
              numchamber+ nummask+
              age2+age3+age4+age5+age6+age7+age8+age9+age10+age11+age12+school2+school3+school4+school5+school6+rankgrade2+rankgrade3+rankgrade4+rankgrade5+rankgrade6+monthsoversea2+monthsoversea3+monthsoversea4+monthsoversea5+monthsoversea6+monthsoversea7+monthsoversea8+monthsoversea9+monthsoversea10+base1+base2+base3+base4+infocenter+honolulu_contact, cw)
summary(chamber1_g, cluster="base")

#ESTIMATE BOOTSTRAPPED SE
boots <- 1000
set.seed(19146)
use.boots1 <- rep(NA, nrow=boots)
use.boots2 <- rep(NA, nrow=boots)
for(b in 1:boots){
  g.star <- cw[sample(1:nrow(cw), replace=TRUE),]
  #First Stage
  boot.first <- lm(use ~ numchamber+ nummask+age2+age3+age4+age5+age6+age7+age8+age9+age10+age11+age12+school2+school3+school4+school5+school6+rankgrade2+rankgrade3+rankgrade4+rankgrade5+rankgrade6+monthsoversea2+monthsoversea3+monthsoversea4+monthsoversea5+monthsoversea6+monthsoversea7+monthsoversea8+monthsoversea9+monthsoversea10+base1+base2+base3+base4+ally_relations+orientation_officers+infocenter+orient_meet+war_interest+honolulu_contact, g.star)
  #Second Stage
  boot.direct <- lm(I(use -
                  coef(chamber1_prepost)['ally_relations']*ally_relations +
                  coef(chamber1_prepost)['orientation_officers']*orientation_officers +
                  coef(chamber1_prepost)['orient_meet']*orient_meet +
                  coef(chamber1_prepost)['war_interest']*war_interest) ~
                  numchamber+ nummask+
                  age2+age3+age4+age5+age6+age7+age8+age9+age10+age11+age12+school2+school3+school4+school5+school6+rankgrade2+rankgrade3+rankgrade4+rankgrade5+rankgrade6+monthsoversea2+monthsoversea3+monthsoversea4+monthsoversea5+monthsoversea6+monthsoversea7+monthsoversea8+monthsoversea9+monthsoversea10+base1+base2+base3+base4+infocenter+honolulu_contact, g.star)
  use.boots1[b] <- coef(boot.direct)["numchamber"]
  use.boots2[b] <- coef(boot.direct)["nummask"]
}

## RESULTS

coef(chamber1_g)["numchamber"]
sd(use.boots1)
coef(chamber1_g)["nummask"]
sd(use.boots2)
nobs(chamber1_g)
AIC(chamber1_g)

coef.vec <- c(-0.026248, -0.02423234)
se.vec <- c(0.0068786, 0.01461281)
var.names <- c("OLS", "Sequential\ng-estimation")

y.axis <- length(var.names):1

layout(matrix(c(2,1),1,2), widths = c(1.5, 5))

par(mar=c(2,2,.5,1))

plot(coef.vec, y.axis, type = "p", axes = F, xlab = "", ylab = "", pch = 19, cex = 2, 
     xlim = c(-.07,.01), xaxs = "r", main = "")
abline(v = 0, col="red", lwd=3, lty=2)

segments(coef.vec-qnorm(.95)*se.vec, y.axis, coef.vec+qnorm(.95)*se.vec, y.axis, lwd =  3)

axis(1, at = seq(-.07,.01,by=.01), labels =  seq(-.07,.01,by=.01), tick = T,
     cex.axis = .8, mgp = c(2,.5,0))
axis(2, at = y.axis, label = var.names, las = 1, tick = T, 
     cex.axis = .8) 

par(mar=c(2,0,.5,0)) 

plot(seq(-0.07,0.01,length=length(var.names)), y.axis, type = "n", axes = F,  xlab = "Estimated Proportion of Respondents Supporting Chemical Weapons Use", ylab = "")

